Mechanism of C-N bonds formation in electrocatalytic urea production revealed by ab initio molecular dynamics simulation

Electrosynthesis of urea from CO2 and NOX provides an exceptional opportunity for human society, given the increasingly available renewable energy. Urea electrosynthesis is challenging. In order to raise the overall electrosynthesis efficiency, the most critical reaction step for such electrosynthesis, C-N coupling, needs to be significantly improved. The C-N coupling can only happen at a narrow potential window, generally in the low overpotential region, and a fundamental understanding of the C-N coupling is needed for further development of this strategy. In this regard, we perform ab initio Molecular Dynamics simulations to reveal the origin of C-N coupling under a small electrode potential window with both the dynamic nature of water as a solvent, and the electrode potentials considered. We explore the key reaction networks for urea formation on Cu(100) surface in neutral electrolytes. Our work shows excellent agreement with experimentally observed selectivity under different potentials on the Cu electrode. We discover that the *NH and *CO are the key precursors for C-N bonds formation at low overpotential, while at high overpotential the C-N coupling occurs between adsorbed *NH and solvated CO. These insights provide vital information for future spectroscopic measurements and enable us to design new electrochemical systems for more value-added chemicals.


Computational Details
DFT Computation Parameters. Further to the computation details outlined in the manuscript, we used the plane-wave cutoff energy of 400 eV and the first order Methfessel-Paxton scheme with a smearing width of 0.2 eV. For static computations, the convergence criteria were 1×10 -5 eV energy differences for solving the electronic wave function while geometry optimization were converged to within 2×10 -2 eV/Å. Interface Modelling. The water/Cu (100) interface was modeled with 32 explicit (four layers) water molecules on a three-layer (3×4×3) Cu(100) surface slabs. The bottom two layers of Cu atoms are fixed at their bulk positions (a=3.67 Å), also water molecules of the topmost layer are kept fixed to keep the water density of the interface close to the density of bulk water. All other atoms are allowed to relax during the simulation.

Molecular Dynamics.
To describe the dynamic nature of hydrogen-bonding network, we carried out ab initio molecular dynamics (AIMD) simulation at 300 K and constructed the free energy profiles. A 1 fs time step with hydrogen mass set to 2 and the convergence criteria for electronic step set to 1×10 -6 eV using the gamma point of the Brillouin zone. This system in canonical ensemble is first heated up from 0K to 300K by rescaling the velocities every 20 steps and then equilibrated at 300K for more than 30ps with a Nose-Hoover thermostat. Then adsorbates were added to the well-equilibrated interface; the initial configurations were determined by minima hopping. 1 The free energy profile for this reaction was obtained from thermodynamic integration of the potential of mean force along the reaction pathway produced either by 2ps AIMD simulations at 12 windows (NO dissociation) or slow-growth with a step size of 0.0008Å for most cases while for hard-to-converge cases, a smaller step size of 0.0005Å or 0.0002Å was used. 2 Both the Eley-Rideal (ER, water as proton source and Langmuir-Hinshelwood (LH, * H at the hollow site as the proton source) reaction mechanisms were considered for the proton coupled electron transfer step. The collective variables (CV) for elementary reactions are defined as the distance between hydrogen and carbon/nitrogen/oxygen atom of the reaction intermediate.

Determination of Electrode Potential
In the present work, we determine the electrode potential range versus SHE by referencing the work function of electrochemical interface to an experimental value of 4.44 eV for the standard hydrogen electrode (SHE).
where is the work function, which can be computed from DFT calculations. Similar schemes were also adopted for the study of electrochemical reduction of CO2. 1, 2, 3, 4 As an example of determining the electrode potential of the electrochemical interface, Figure S2 demonstrates that the work function varies with time along the AIMD trajectories.
We found that although the temperature and potential energy converge quickly ( Figure S1), it takes about 24ps for work functions to reach an acceptable equilibrium. The work function is reduced from the initial configuration until it reaches the equilibrium. However, computing work functions for each point of the trajectories requires significant computational resources.
For other cases, we could not sample all the points and instead, only dozens of the last few picoseconds were considered. The work function values for Cu (100) with various intermediates are summarized in Table S1 and Table S2. These data show that although the

Constant Potential Corrections
All the calculations were done when the number of electrons was fixed, which meads the work function as well as the electrode potential referenced with standard hydrogen electrode changes along the reaction coordinates. However, the electrochemical measurements were conducted under fixed potential referenced with a certain reference electrode. To get electrochemical barrier at constant potential, we adopted methods developed by Chan and Nørskov based on a capacitor model: q  , and   corresponds to the energy correction due to change of electrode potential, charges and workfunctions. 5,6 In this present work, we calculated the capacitors instead according to the following equation 3, which also proved valid. 4 Using the method proposed in Ref. 13, the calculated capacitance (C) is 1.27 e/V. When the electrode potential values in Table S1 and S2 are assigned to the parameters in eq.3, we found that for most cases (only one exception), the energy correction term is less than 0.06 eV (corresponding to the work function change of 0.3eV), which is less than the standard deviation of the calculated barriers in Table S3-8. Hence, we concluded that the constant potential correction terms are minor in the present study.